! SK INLET/OUTLET 1D CHARACTERISTIC NONREFLECTING BCS  



subroutine usdef1( un, u,  ub, ucb, ucsb  &

                     , pn,p, pb,pcb, pcsb &

                     , ron,rho,rhob,rhocb,rhocsb &

                     ,ds1,esm,dt,dtmin,epsilon2,p00,tout,ubg1,finlet1)

	!IMPLICIT REAL *8(a-h,o-z)

use Data_cfoam, only: epsMaximumPrinciple

implicit none



    real *8 esm(3)

	real *8  un(3), u(3),  ub(3), uf(3), ucb(3), ucsb(3), ucf(3), ucsf(3)

	real *8 pn,p, pb, pf, pcb, pcsb, pcf, pcsf

	real *8 rhon,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf

	real *8 ds1,ds2,dt,sound,rofon,dtmin,tout

	real *8 ut(3,9)

! SK EXTRA PARAMETERS SPECIAL FOR THE OUTLET

    real *8 uf1,uf2,uf3,pout,finlet1



! SK EXTRA VARIABLES

    real *8 qb,qcb,qcsb,rf,rcf,rcsf,rwellb

	real *8 qwellb,qmaxb,qminb,rwellf,rmaxf,rminf

	real *8 gl,gr,rnr,rnl,qnl,qnr,char1f,char1b

	real *8 smachl,smachr,smach,gn,ubg1(3),ubg(3)

! SK



! END EXTRA PARMS

    real *8 ei,z1,vtmp,dttmp

    real *8 char3,q,r,qn,rn,u1x,ron,vf,wf,vb,wb

    real *8 char3f,soundf,gf,sf,qf,qcf,qcsf,char2f,qwell,qmaxf,qminf

    real *8 char3cf

    real *8 char3csf

    real *8 char3b,soundb,gb,sb,rb,rcb,rcsb,char2b,rwell,rmaxb,rminb

    real *8 char3cb

    real *8 char3csb

    real *8 charc1,charc2,charc3

    real *8 a(3,3)

    real *8 epsilon

    real *8 es(3)

    real *8 defb

    real *8 alfp,dtmint,gam,gam2,gamx,gamy,p00,roout,s,scb,scsb,sn,cv

    integer j

    

    external defb

    real *8 epsilon2



    real(8) :: sn1,sn2, ut12,ut13,ut22,ut23, pn1,pn2, ron1,ron2

    real(8) :: alfa



    epsilon=epsilon2

    !if (pcsb.lt.0.95d5) epsilon=0.2

    !if (pcsb.gt.1.15d5) epsilon=0.2

    gam=1.4d0

    gamx=2.0d0*dsqrt(gam)/(gam-1.0d0)

    gamy=((gam-1.0d0)/(2.0d0*gam))

    gam2=2.0d0*gam

    es=esm/dsqrt(esm(1)**2+esm(2)**2+esm(3)**2)

!	p00=101325.0d0 

	cv=717.65d0





	call coor(es(1),es(2),es(3),a)

! 1-u 2-ub 3-uf 4-ucb 5-ucsb 6-ucf 7-ucsf

    call  transform(u,ut(1,1),a)

    call  transform(ub,ut(1,2),a)





    call  transform(ucb,ut(1,4),a)

    call  transform(ucsb,ut(1,5),a)



!    ubg(1)=ubg1

!    ubg(2)=0.0

!    ubg(3)=0.0



    ubg(1:3)=ubg1(1:3)



! SK TRANSFORM THE BACKGROUND VELOCITY

    call  transform(ubg,ut(1,3),a)

!if(finlet1.eq.1.0d0)then

!    ut(1:3,3) = -ubg(1:3)

!end if

! END OF TRANSFORM BCG VELOCITY 



    char3 =ut(1,1)

    char3b=ut(1,2)

 

    char3cb=ut(1,4) 

    char3csb=ut(1,5)





! SK SWITCHES BETWEEN OUTLET/INLET -- SAME AS IN THE PREVIOUS VERSION////////////////////////

	alfp=0.0

!	alfp=1.0

!    if(char3csb+ut(1,3).lt.0)alfp=1.0d0

!    call GetAlfa(char3csb+ut(1,3),alfp)

	if(finlet1.eq.1.0d0) alfp=1.0d0

	pout=p00 !*alfp+pcsb*(1.0d0-alfp) ! outside pressure

	roout=p00/(tout*cv*(Gam-1.0d0))

	

!    if(ut(1,3).gt.0)alfp=0.0d0



! SK //////////////////////////////////////////////



!    char3 =ut(1,1)*(1.0d0-alfp) + ut(1,3)*alfp

!    char3b = ut(1,2)*(1.0d0-alfp) + ut(1,3)*alfp

!    char3cb = ut(1,4)*(1.0d0-alfp) + ut(1,3)*alfp

!    char3csb = ut(1,5)*(1.0d0-alfp) + ut(1,3)*alfp





    s=dlog(p)-gam*dlog(rho)

    sb=dlog(pb)-gam*dlog(rhob)



    scb=dlog(pcb)-gam*dlog(rhocb)

    scsb=dlog(pcsb)-gam*dlog(rhocsb)



    gb=gamx*dexp(scsb/gam2)

    soundb=dsqrt(gam*pcsb/rhocsb)



! SK COMPUTE THE REFERENCE FORWARD STATE



   uf1=ut(1,5)*(1.0d0-alfp)+ut(1,3)*alfp

   uf2=ut(2,5)*(1.0d0-alfp)+ut(2,3)*alfp

   uf3=ut(3,5)*(1.0d0-alfp)+ut(3,3)*alfp

   rhof=rhocsb*(1.0d0-alfp)+roout*alfp

   pf=pcsb*(1.0d0-alfp)+pout*alfp



!   uf1=ut(1,3) - (1.0d0-alfp) * 2d0 * ut(1,3)

!   uf2=ut(2,3) - (1.0d0-alfp) * 2d0 * ut(2,3)

!   uf3=ut(3,3) - (1.0d0-alfp) * 2d0 * ut(3,3)

!   rhof= roout

!   pf= pout



   sf=dlog(pf)-gam*dlog(rhof)

   gf=gamx*dexp(sf/gam2)

   soundf=dsqrt(gam*pf/rhof)

   char3f=uf1





! END FORWARD STATE





! SK THEN JUST COPY THE RELEVANT PART FROM THE MAIN USDEF ROUTINE (WITH ALL EPS=THE SAME)



    qf=char3f-gf*pf**gamy

    qnl=qf



!SK  + 4 TERMS IN THE OTHER DIRECTION  (B)

    q=char3-gb*p**gamy      

    qb=char3b-gb*pb**gamy

    qcb=char3cb-gb*pcb**gamy

    qcsb=char3csb-gb*pcsb**gamy

! SK + EXTRAPOLATION AND CORRECTION IN THE OTHER DIRECTION (B)



    char2b=char3csb-soundb	

	qwellb=(qcsb-qcb)/(0.5d0*dt)+char2b*(q-qb)/ds1

	qmaxb=dmax1(q,qcb,qb)+dt*qwellb + epsMaximumPrinciple

	qminb=dmin1(q,qcb,qb)+dt*qwellb - epsMaximumPrinciple

	qnr=(2d0*qcsb-qb*(1.0d0-epsilon))/(1.0d0+epsilon)

	if(qnr>qmaxb) qnr=qmaxb

	if(qnr<qminb) qnr=qminb



! SK



	r=char3+gb*p**gamy

	rb=char3b+gb*pb**gamy

	rcb=char3cb+gb*pcb**gamy

	rcsb=char3csb+gb*pcsb**gamy



	char1b=char3csb+soundb

	rwellb=(rcsb-rcb)/(0.5d0*dt)+char1b*(r-rb)/ds1	

	rmaxb=dmax1(r,rcb,rb)+dt*rwellb + epsMaximumPrinciple

	rminb=dmin1(r,rcb,rb)+dt*rwellb - epsMaximumPrinciple



    rnr=(2d0*rcsb-rb*(1.0d0-epsilon))/(1.0d0+epsilon)

	if(rnr>rmaxb) rnr=rmaxb

	if(rnr<rminb) rnr=rminb

!SK +4 TERMS IN THE OTHER DIRECTION (F)

        rf=char3f+gf*pf**gamy

        rnl=rf





!  SK INVARIANTS DEFINITION END





! SK ALL CHAR'S AT THE RHS BELOW MUST BE REFERRED TO THE CENTRES 

!    sound=soundf+soundb        ! not used VSemiletov

!    charc3=char3f+char3csb

!    charc1=char3+sound   

!    charc2=char3-sound         ! not use, VSemiletov







! SK DEFINITION //////////////////////////////////////

    gr=gb

    gl=gf

!    char3f=char3csf  ! not use, VSemiletov

!    char3b=char3csb  ! not use, VSemiletov

! END DEFINITION//////////////////////////////////////



! SK THE BASE CASE: SUBSONIC LEFT/RIGHT GOING WAVES//////////////////////



    gf=gl

    gb=gr

    rn=rnr

    qn=qnl



    u1x=qn+gf*(rn-qn)/(gf+gb)

!    u1x = u1x*(1.0d0-alfp) + ut(1,3)*alfp

    call GetAlfa(u1x, alfa)

    if(alfa>=1.0d0)then

        sn=sf

        ut(1,8)=0.0d0

        ut(2,8)=uf2

        ut(3,8)=uf3

        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))

        ron=(pn/exp(sn))**(1.0d0/gam)

    else if(alfa<=0.0d0)then

        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)

        ut(1,8)=0.0d0

        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)

        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)

        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))

        ron=(pn/exp(sn))**(1d0/gam)

    else

        sn1=sf

        ut12=uf2

        ut13=uf3

        pn1=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))

        ron1=(pn1/exp(sn1))**(1.0d0/gam)

        

        sn2=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)

        ut22=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)

        ut23=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)

        pn2=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))

        ron2=(pn2/exp(sn2))**(1d0/gam)

        

        sn = alfa * sn1 + (1.0d0-alfa) * sn2

        ut(1,8)=0.0d0

        ut(2,8) = alfa * ut12 + (1.0d0-alfa) * ut22

        ut(3,8) = alfa * ut13 + (1.0d0-alfa) * ut23

        pn = alfa * pn1 + (1.0d0-alfa) * pn2

        ron = alfa * ron1 + (1.0d0-alfa) * ron2

    end if

    

!    if(u1x<=0.0) then   ! if(charc3<=0.0) then ! SK changed, 2011-05-17

!        sn=sf

!        ut(1,8)=0.0d0

!        ut(2,8)=uf2

!        ut(3,8)=uf3

!        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))

!        ron=(pn/exp(sn))**(1.0d0/gam)

!        !if (k.eq.5051) write(14,*) 'f',pn,ron

!        !if (k.eq.5051) write(14,*) 'fqr',q,r

!    else 

!        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)

!        ut(1,8)=0.0d0

!        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)

!        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)

!        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))

!        ron=(pn/exp(sn))**(1d0/gam)

!        !if (k.eq.5051) write(14,*) 'b',pn,ron

!        !if (k.eq.5051) write(14,*) 'bqr',q,r

!    end if



! SK END OF SUBSONIC ////////////////////////////////////////////



! SK RIGHT-GOING SUPERSONIC/////////////////////////////////////////////





    if(char1b>0.0 .and. char2b>0.0) then

        gb=gr

        gf=gb

        rn=rnr

        qn=qnr



        u1x=0.5d0*(rn+qn)

        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)

        ut(1,8)=0.0

        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)

        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)

        pn=((rn-qn)/2.0d0/gb)**(2.0d0*gam/(gam-1.0d0))

        ron=(pn/exp(sn))**(1d0/gam)

    end if





! SK END OF R/G SUPERSONIC ////////////////////////////////



! SK LEFT-GOING SUPERSONIC/////////////////////////////////////////////



    if(char1b<=0.0 .and. char2b<0.0) then

        gb=gl

        gf=gb

        rn=rnl

        qn=qnl



        u1x=0.5d0*(rn+qn)

        sn=sf

        ut(1,8)=0.0

        ut(2,8)=uf2

        ut(3,8)=uf3

        pn=((rn-qn)/2.0d0/gf)**(2.0d0*gam/(gam-1.0d0))

        ron=(pn/exp(sn))**(1.0d0/gam)

    end if



! SK END OF L/G SUPERSONIC ////////////////////////////////

   



    call transformr(ut(1:3,8),ut(1:3,9),a)

	do j=1,3    

        un(j)=ut(j,9)+u1x*es(j)

	end do



	vtmp=dsqrt(un(1)**2+un(2)**2+un(3)**2)

    sound=dsqrt(gam*pn/ron)

    dtmint=ds1/(sound+vtmp)

	if (dtmint.lt.dtmin) dtmin=dtmint ! definition of Time step



    return      

end subroutine usdef1